##############################
# STATE REACH AND DEVELOPMENT // PSRM // Replication
#
# Plotting the input data to 
# the state reach measure
#
##############################

# African cowcodes
afro.cow <- c(400:626, 651)
afro.cow <- afro.cow[!afro.cow %in% c(581, 402, 403, 590, 591)] ## Only mainland Africa!!


# Data Overview: Admin Unit Maps ####


# Load data
admin.shp <- readOGR(dsn = file.path("data", "africa_regions_panel.GeoJSON"), 
                     layer="africa_regions_panel")

# FIGURE A2: Map Admin units in 1 point in time

# ... set year of the map
map.yr <- 2016

# ... select this map
this.map <- admin.shp[admin.shp$startyr <= map.yr &
                        admin.shp$endyr >= map.yr,]

# ... combine to countries for thick borders
cow.map <- gUnaryUnion(this.map, id = this.map$cowcode)
cow.map <- remove.holes(cow.map, max.size = 1e6)

# ... plot
png(paste0(fig.path, "figa2_data_admunits_map.png"), width = 7, height = 7, units = "in", res = 300)
par(mar = c(0,0,0,0))
plot(this.map, border = "red")
plot(cow.map, border = "black", lwd = 2, add = T)
dev.off()

# FIGURE A1: Plot Number of units in data set

# ... to yearly obs
year.df <- unlist(lapply(c(1:length(admin.shp)), function(i){
  admin.shp$startyr[i]:admin.shp$endyr[i]
}))
year.df <- data.frame(year = year.df, obs = 1)
year.df <- aggregate.data.frame(list(obs = year.df$obs),
                                list(year = year.df$year), FUN = sum)

# ... plot
png(paste0(fig.path, "figa1_data_admunits_obs.png"), width = 3, height = 2.5, units = "in", res = 300)
g <- ggplot(year.df, aes(x = year, y = obs)) + geom_line() +
  xlab("Year") + ylab("Number of regions") +
  theme_minimal()
print(g)
dev.off()


# FIGURE A1: Plot observed changes by year

# ... get first year by country
first.yr.df <- lapply(c(1:length(admin.shp)), function(i){
  admin.shp$startyr[i]:admin.shp$endyr[i]
})
first.yr.df <- data.frame(cowcode = rep(admin.shp$cowcode, unlist(lapply(first.yr.df, length))),
                          year = unlist(first.yr.df))
first.yr.df <- aggregate.data.frame(list(first.yr = first.yr.df$year),
                                list(cowcode = first.yr.df$cowcode), FUN = min)

# .. make change df
change.df <- join(admin.shp@data, first.yr.df)
change.df <- change.df[change.df$startyr != change.df$first.yr,]

# ... count nubmer of changes per year
change.df$change.num <- 1
change.df <- aggregate.data.frame(list(obs = change.df$change.num),
                                  list(year = change.df$startyr), FUN = sum)

# ... plot
png(paste0(fig.path, "figa1_data_admunits_changes.png"), width = 3, height = 2.5, units = "in", res = 300)
g <- ggplot(change.df, aes(x = year, y = obs)) + geom_line() +
  xlab("Year") + ylab("Number of changes") +
  theme_minimal()
print(g)
dev.off()


# FIGURE A2: Plot observed number by country

# ... make data frame
all.yr.df <- lapply(c(1:length(admin.shp)), function(i){
  admin.shp$startyr[i]:admin.shp$endyr[i]
})
all.yr.df <- data.frame(cowcode = rep(admin.shp$cowcode, unlist(lapply(all.yr.df, length))),
                        countryiso = rep(admin.shp$countryiso, unlist(lapply(all.yr.df, length))),
                        year = unlist(all.yr.df),
                        obs = 1)
all.yr.df <- aggregate.data.frame(list(obs = all.yr.df$obs),
                                    list(cowcode = all.yr.df$cowcode, 
                                         countryiso = all.yr.df$countryiso,
                                         year = all.yr.df$year), FUN = sum)

# ... plot
png(paste0(fig.path, "figa2_data_admunits_obsbycow.png"), width = 7, height = 12, units = "in", res = 300)
g <- ggplot(all.yr.df, aes(x = year, y = obs)) + geom_line() + 
  facet_wrap(~countryiso, ncol = 4) + xlab("Year") + ylab("Number of regions") +
  theme_minimal()
print(g)
dev.off()

# FIGURE A3: Data Overview: Michelin Maps ####


# ... load Michelin map extents
ext.sp <- readOGR(dsn = file.path("data", "michelin_extents.GeoJSON"), 
                  layer="michelin_extents")

# ... cshapes 2016
cshp.2016 <- cshp(date = as.Date("2016-01-01"))
cshp.2016 <- cshp.2016[cshp.2016$gwcode %in% afro.cow,]

# ... colors
col.fill = col2rgb(c("red","blue","green"), alpha = .5) / 255

# ... plot
png(paste0(fig.path, "figa3_data_michelin_spatcov.png"), width = 5, height = 5, units = "in", res = 300)
par(mar = c(0,0,0,0))
plot(as_Spatial(cshp.2016))

for(i in seq_along(unique(ext.sp$type))){
  plot(ext.sp[ext.sp$type == unique(ext.sp$type)[i],], 
       col = rgb(col.fill[1,i],col.fill[2,i],col.fill[3,i],.5), add = T)
}
dev.off()

# FIGURE A4: Temporal Coverage 

# ... load data
year.df <- read.csv(file.path("data", "michelin_tempcov.csv"),
                    stringsAsFactors = F)

# ... map to factor
year.df$map.fac <- as.numeric(factor(year.df$map, levels = unique(year.df$map)))

# ... plot
png(paste0(fig.path, "figa4_data_michelin_tempcov.png"), width = 5, height = 2.5, units = "in", res = 300)
g <- ggplot(year.df, aes(y = map.fac, x = year)) + geom_point(pch = 15) +
  geom_hline(yintercept = 1.5, lty = 2) + 
  scale_y_continuous(name = NULL, 
                     breaks = unique(year.df$map.fac), 
                     labels = unique(year.df$map)) +
  scale_x_continuous(name = "Year", 
                     breaks = c(1966, seq(1970, 2010, 10), 2017))  +
  theme_minimal()
print(g)
dev.off()

# Cleanup large files
rm(list = c("admin.shp"))
